/*
    Copyright 2007-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    Definitions and explicit instantiations of the ir_grid class
    members. */

// $Id$

#include <CCfits/CCfits>
#include "ir_grid.h"
#include "grid.h"
#include "arepo_grid.h"
#include "arepo_emission.h"
#include "blitz-fits.h"

#ifdef WITH_AREPO
extern "C" {
#include "allvars.h"
}
#endif

using namespace blitz;

/** Constructor loads grid structure and cell data from a FITS file.
    The grid structure is loaded from a GRIDSTRUCTURE HDU by the base
    class constructor, and then the auxiliary quantities are loaded
    from the data hdu. */
template <template<typename> class T_grid, 
	  template<typename> class sampling_policy, 
	  typename rng_policy>
mcrx::ir_grid<T_grid, sampling_policy, rng_policy>::
ir_grid (boost::shared_ptr<T_grid_impl> g,
	 CCfits::ExtHDU& structure_hdu,
	 CCfits::ExtHDU& data_hdu,
	 const density_generator& dg,
	 int n_lambda) :
  emission_collection<polychromatic_policy, sampling_policy, rng_policy>(),
  g_(g), n_lambda_(n_lambda)
{
  structure_hdu.readKey("lengthunit", units_ ["length"]);

  // get the amount of dust in cells.  For this we need the density generator
  const int nc = this->n_cells();
  m_dust_.resize(nc, dg.make_density_vector(0,0).size());

  // read amount of metals in cells
  array_1 m_gas, m_metal;
  read (data_hdu.column("mass_gas"), m_gas);
  read (data_hdu.column("mass_metals"), m_metal);

  // loop over cells and create dust mass array and emitter objects
  int i = 0;
  for (typename T_grid_impl::iterator c = this->begin (), e = this->end(); 
       c != e; ++c, ++i) {
    // density generator takes and returns DENSITIES, but we want MASS
    m_dust_(i, Range::all()) = 
      dg.make_density_vector(m_gas(i)/c->volume(), m_metal(i)/c->volume()) *
      c->volume();
  }

  make_emitters();
}

// helper function to call the arepo_emission method if that's the type
template<typename T> void update_bounding_boxes(const T&) {assert(0);};
template<typename T> void 
update_bounding_boxes(const mcrx::arepo_grid<T>& g) {
  T::update_bounding_boxes(g); };


/** Constructor getting the data from the Arepo data structures. */
template <template<typename> class T_grid, 
	  template<typename> class sampling_policy, 
	  typename rng_policy>
mcrx::ir_grid<T_grid, sampling_policy, rng_policy>::
ir_grid (boost::shared_ptr<T_grid_impl> g,
	 const density_generator& dg,
	 int n_lambda) :
  emission_collection<polychromatic_policy, sampling_policy, rng_policy>(),
  g_(g), n_lambda_(n_lambda)
{
#ifdef WITH_AREPO

  // First copy units from the arepo units
  units_["length"]=arepo::arepo_units["length"];
  units_["time"]=arepo::arepo_units["time"];
  units_["mass"]=arepo::arepo_units["mass"];

  // get the amount of dust in cells.  For this we need the density generator
  const int nc = this->n_cells();
  m_dust_.resize(nc, dg.make_density_vector(0,0).size());

  const T_float convert_density=
    All.UnitDensity_in_cgs*units::convert("g/cm^3",units_["mass"]+"/"+units_["length"]+"^3");

  int i = 0;
  for (typename T_grid_impl::iterator c = this->begin (), e = this->end(); 
       c != e; ++c, ++i) {

    assert(c->volume()>0);
    const T_float rho_gas = SphP[i].Density*convert_density;
    const T_float rho_met = rho_gas*P[i].Metallicity;

    // density generator takes and returns DENSITIES, but we want MASS
    m_dust_(i, Range::all()) = 
      dg.make_density_vector(rho_gas, rho_met) * c->volume();
  }

  make_emitters();
  // since this is the arepo grid, we also now need to set up the
  // bounding boxes of the emitters
  update_bounding_boxes(*g_);

#else
  cerr << "Not compiled with Arepo support" << endl;
  exit(1);
#endif
}

/** The ir_grid is not like a conventional grid in that the quantities
    are stored in the grid cells. Because the main function of this
    grid is to calculate dust emission, those quantities are stored in
    the grid object as an array, not in the cells. This means that we
    are unable to build this grid the traditional way using a grid
    factory. Instead this constructor uses a structure vector and a
    mass vector. */
template <template<typename> class T_grid,
	  template<typename> class sampling_policy, 
	  typename rng_policy>
mcrx::ir_grid<T_grid, sampling_policy, rng_policy>::
ir_grid (boost::shared_ptr<T_grid_impl> g,
	 const array_2& mass,
	 const T_content& lambda,
	 const T_unit_map& units) :
  emission_collection<polychromatic_policy, sampling_policy, rng_policy>(),
  g_(g), m_dust_(mass), n_lambda_(lambda.size())
{
  units_["length"] = units.get("length");
  // Because the dust mass is already supplied, we just need to set up
  // the emitters.
  make_emitters();
}


template <template<typename> class T_grid, 
	  template<typename> class sampling_policy, 
	  typename rng_policy>
void
mcrx::ir_grid<T_grid, sampling_policy, rng_policy>::
make_emitters()
{
  const int nc = this->n_cells();

  // Allocate SED array
  cout << "Allocating a (" << nc << ", " << n_lambda_
       << ") memory block for cell emission, "
       << 1.0*nc*n_lambda_*sizeof(T_float)/(1024*1024*1024) << " GB.\n";
  L_lambda_.resize(nc, n_lambda_ );
  // Set to invalid value
  L_lambda_ = blitz::quiet_NaN(T_float());

  int i = 0;
  array_1 temp_lum;
  for (typename T_grid_impl::iterator c = this->begin (), e = this->end(); 
       c != e; ++c, ++i) {

    // create an emitter object referencing the main L_lambda array
    temp_lum.weakReference(L_lambda_(i, Range::all ()));
    assert (temp_lum.isStorageContiguous());
    
    typename T_grid_impl::T_data* dd = new T_emitter (temp_lum, *c);
    c->set_data(dd);
  }
}


// Explicit instantiations for the two cases we use. If someone starts
// using rejection sampling or invents a new grid, these need to be
// added here.
namespace mcrx {
  template class ir_grid<adaptive_grid, cumulative_sampling, local_random>;
#ifdef WITH_AREPO
  template class ir_grid<arepo_grid, cumulative_sampling, local_random>;
#endif
}
